The low temperature Fulde-Ferrell-Larkin-Ovchinnikov phases in 3 dimensions 
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We consider the nature of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phases in three dimen- 
sions at low temperature. We introduce a new method to handle the quasiclassical equations for 
superconductors with space dependent order parameter, which makes use of a Fourier expansion. 
This allows us to show that, at T = 0, an order parameter given by the linear combination of three 
cosines oscillating in orthogonal directions is preferred over the standard single cosine solution. The 
transition from the normal state to this phase is first order, and quite generally the transition below 
the tricritical point to the FFLO phases is always first order. 

PACS numbers : 74.20.Fg, 74.25.Op, 74.70.Kn, 74.81.-g 
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Although they have been introduced quite a long 
time ago [1,2], Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
phases are still of very high interest. On the experi- 
mental side they should appear as the superconducting 
phases in extremely high field superconductors, which 
are obviously of very high practical interest. On the the- 
oretical side, quite surprisingly, a precise understanding 
of the order paramater and nature of the transitions in 
these phases has not yet been reached, and this question 
is naturally of essential importance for the experimen- 
tal identification of these phases, which have not been 
observed to date unambiguously. These FFLO phases 
correspond to a spontaneous symmetry breaking of the 
standard BCS superfluid phase in the presence of an ef- 
fective field, which has the only effect of inducing a differ- 
ence in chemical potential between the two populations 
involved in the formation of Cooper pairs. This leads to 
an inhomogeneous superfluid with a space dependent or- 
der parameter. The chemical potential difference is due 
for example to the coupling of the effective field with the 
spins of the particles. Naturally in standard supercon- 
ductors a real magnetic field gives rise to a coupling to 
the electronic orbital degrees of freedom which is much 
stronger than the coupling to the spins and is ordinarily 
responsible for the critical field. Here we investigate the 
original problem considered by Larkin and Ovchinnikov 
where this orbital coupling is absent or essentially negli- 
gible. Accordingly we consider the simplest model with 
a spherical Fermi surface. 

In their original work Larkin and Ovchinnikov [2] (LO) 
explored at T = the possibility of a second order phase 
transition and looked, in the three-dimensional case, for 
the order parameter corresponding to the lowest energy. 
They found that the stablest state has a simple one- 
dimensional space dependence A(r) ~ cos(q.r) for the 
order parameter. Nevertheless they left open the possi- 
bility that a first order transition could exist at higher 
fields. However looking for a first order transition is a 
much more complicated problem than exploring a sec- 
ond order one. Indeed LO could solve for the second 



order transition by performing an expansion of the free 
energy up to fourth order in powers of the order parame- 
ter, while no expansion is in principle allowed in the case 
of a first order transition. Similarly LO could restrict the 
space dependence to a sum of plane waves, while for first 
order this is no longer correct. In this paper we consider 
the original problem raised by LO and show, by solving 
the quasiclassical equations with the use of a Fourier ex- 
pansion, that the transition is always first order. We find 
in particular that, at T = 0, the preferred order param- 
eter has a cubic symmetry and is surprisingly very near 
one of the order parameters investigated explicitely by 
LO. 

The question of the order of the transition to the FFLO 
phases has given rise recently to a good deal of work. 
Burkhardt and Rainer [3] used quasiclassical equations 
to show that, for the two-dimensional case, the transi- 
tion to A(r) ~ cos(q.r) stays second order all the way 
up to the tricritical point, where the FFLO transition 
line meets the standard BCS line. Very recently we have 
investigated [4] the specific nature of the resulting or- 
der parameter in the low temperature regime. On the 
other hand in the three- dimensional case, Matsuo [5] et 
al showed, again by solving quasiclassical equations, that 
the transition to this same order parameter becomes first 
order at low temperature, above T/T c q = 0.075, where 
T c q is the standard critical temperature for equal spin 
chemical potentials. It stays then first order up to the 
tricritical point. A different line of attack has been to 
investigate this problem in the vicinity of this tricriti- 
cal point (located at T tcv /T C Q = 0.561). Indeed because 
of the proximity of the second order transition to the 
standard BCS phase with uniform order parameter, an 
expansion up to sixth order [6-8] in the order parameter 
can always be used. This has been done by Houzet et al 
[7] who examined a number of possible order parameters 
and concluded that the best one was indeed the simple 
cosine A(r) ~ cos(q.r). Quite recently we have inves- 
tigated this same problem analytically [8] and came to 
the same conclusion by showing that the states with the 
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smallest number of plane waves are favored. 

These explorations brought a somewhat surprising in- 
formation. In the vicinity of the tricritical point [7,8] the 
location of the first order transition is extremely close 
to the standard FFLO second order transition. This re- 
mains true [5], for the transition to the simple cosine, 
down to the lowest temperatures. This feature suggests 
that the first order transition is not very far from being 
second order. If this were the case, this would allow an 
expansion in powers of the order parameter (for example 
up to sixth order) to remain valid at any temperature, 
and not only near the tricritical point. This has been 
explored by Houzct et al [7], but this approach does not 
work because somewhat below the tricritical point the 
coefficient of the sixth order term changes sign. Actually 
one can show that, whatever the order of the expansion, 
the coefficients are changing sign when the temperature 
is lowered. So an expansion in powers of the order pa- 
rameter can not succeed. This leads to the conclusion 
that the quasiclassical approach is necessary to deal with 
this problem. On the other hand another feature valid 
near the tricritical point can be used. Indeed our explo- 
ration [8] showed that the order parameter at the first 
order transition is very near the simple plane waves su- 
perposition coming in the LO analysis of the second order 
phase transition. Hence we may expect that a Fourier ex- 
pansion will be a good approach at any temperature. In 
this paper we show indeed that introducing a Fourier ex- 
pansion in the quasiclassical equations is a very efficient 
scheme to deal with this non-linear problem. We will see 
that such an expansion converges very rapidly toward 
the exact result and that, as a consequence, the first few 
terms give an excellent approximation. 

Here we will use Eilenberger's original equations which 
are enough to handle our problem. However an essen- 
tial strength of the quasiclassical method is its ability 
to be generalized for full inclusion of many body effects 
[9]. The quasiclassical Green's functions g(u),\s., r) = 
— J d£kG(u>, k, r) are obtained by " ^-integrating " the 
usual temperature Green's functions G(w„,k, r). Here 
£k is the kinetic energy measured from the average Fermi 
level (1/2) (/Xf + MJ.)- The frequency ui will turn at some 
stage to be the Matsubara frequency ui n — irT(2n + 1). 
The existence of the FFLO phases is due to the effective 
field p,= {fi-\ — /Uj.)/2 produced by the difference between 
spin up and spin down chemical potentials. It will merely 
appear [2] by the replacement of u> n by u> n = u„ — ip, in 
the calculations. We simplify the writing by taking K = 1 
and m — 1/2. Eilenberger's equations for the off-diagonal 
f(u>, k, r) quasiclassical propagators then read [10] : 

(w + k.V)/(w,k,r) = A(r) 5 (w,k,r) 
(w - k.V)/+(w, k, r) = A*(r) fl (w, k, r) (1) 

The diagonal propagator g is given in terms of / and f + 
by the normalization condition : 



g(w, k, r) = (1 - f(w, k, r)/ + (w, k, v)) 1 ' 2 (2) 

from which one can [10] deduce easily an equation for g 
similar to Eq.(l). Actually the equations Eq.(l) for / and 
/+ are related since [10] /*(— w,k,r) = f*(w, — k, r) = 
f + (uJ 7 k, r) with analogous relations holding for g(u), k, r). 

We will consider here the very wide class of order pa- 
rameters which admit a three-dimensional Fourier expan- 
sion. The order parameters considered by LO fall in par- 
ticular in this class. We will also assume that the order 
parameter is real. It has indeed been shown [8] that this 
class is favored near the tricritical point, and the order 
parameters considered specifically by LO are also real. 
For clarity we will only consider explicitely [11] the case 
of a simple cosine order parameter A(x) = 2Acos(<;x) 
(the value of q and A are ultimately obtained by mini- 
mization, see below) and then explain the modifications 
occuring when one goes to the general case. For fixed k 
Eilenberger's equations are a set of first order equations 
for the variation of the Green's functions along k. So 
we take a reduced variable along this direction by set- 
ting r = kA which gives A(x) = 2Acos(QA) where we 
have introduced Q = kpq cos 6 with 9 the angle between 
k and the x axis. Then we make a Fourier expansion of 
the Green's functions : 

f(X) = J2f n e inQX (3) 

n 

and similarly for / + (A) and g(X). Parity implies 
f(—X) = f+(X) and g{-X) = g(X), which gives 
f* n = f n and g- n = g n . Inserting these expansions 
in Eq.(l) we find the recursion relations: 

nQA 

dn --^ + (nQy {9n - 1+9n+l) 
A 

9n = — ~{d n -i +d n+1 ) (4) 
nQ 

where we have set d n — (/„ — f£)/2i, from which 
fn = (i — oj/nQ)d n is obtained. These equations show 
that 52 P +i = and d 2p = 0. Eq.(4) are linear and must 
be supplemented by the normalization condition Eq.(2). 
The n — component is enough and it provides us pre- 
cisely with the spatial integral go = J drg(iv 7 k, r) needed 
to calculate the free energy : 

0O = l-X>»S-n + /™/±J ^ 

The critical temperature T is then obtained by writing 
the equality Q s — Q n between normal and superfluid 
state free energies. When this is done from Eilenberger's 
expression [12] for Cl s , we find for our simple cosine, af- 
ter some transformations [11], that ln[T/T sp (p,/T)] is the 
minimum of: 
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~ "aX / dujRc l / dug (u>, qu)-l + — ] (6) 

with a; = lu — ip and where T sp (p/T) is the location of 
the spinodal transition, which is the simple continuation 
of the standard BCS transition. Since for homogeneity 
the above quantity is only a function of T/p, A/p and 
q/p one has to minimize with respect to A/p and q/p, at 
fixed T/p. In this way one finds the critical temperature 
T as a function of T/p, hence T as a function of p. 

Now the interesting point is the large n behaviour of 
g n and d n . The general solution of the recursion rela- 
tions Eq.(4) is a linear combination of two independent 
solutions. For large n the recursion relation for g n simpli- 
fies into A 2 (g n+2 + g n -2) + (nQ) 2 g n = 0. This equation 
has very rapidly growing solutions behaving as g 2 p+2 ~ 
(-1) P (2Q/A) 2p (p!) 2 . Naturally these solutions arc not 
physically acceptable. On the other hand the recursion 
relation has also a solution satisfying g n +2 <C g n <C g n -i 
and behaving as g 2p ~ (-l) p (A/2Q) 2 J>(l/p!) 2 , which is 
the physical solution we are looking for. This solution is 
found only if g and g 2 are related by a specific boundary 
condition. 

Up to now no approximation has been made. Now the 
very fast decrease of g n and d n provides an easy way to 
obtain a set of approximate solutions, which moreover 
converges rapidly to the exact one, all the more since 
these are <? 2 and c? 2 which come in Eq.(5) for the calcu- 
lation of g . Since g n and d n are very small for large 
n we just take g n and e?„ +1 to be zero for some fixed 
value n = N + 1 and beyond. This serves as boundary 
condition. Then we work backward to obtain the whole 
set of Fourier components and normalize them properly 
through the normalization condition Eq.(5). When we 
let N — > oo we find the exact result for go. The recursion 
relations Eq.(4) are very convenient and very fast for a 
numerical implementation and in practice the situation 
is not very different from having an analytical expression 
for g n . Basically the solution of Eq.(l) is reduced to some 
simple algebra. The simplest of the set of converging ap- 
proximations corresponds to take N = 1 and it is given 
cxplicitely by : 



go = [1 + 2A 2 



(7) 



(w 2 + Q 2 ) 2J 

This is already a quite non trivial approximation. Since 
it is correct up to order A 2 it gives the proper location for 
the standard FFLO second order transition line. More- 
over it gives qualitatively and semiquantitatively the cor- 
rect results, with a first order transition down from the 
TCP which becomes a second order transition at low tem- 
perature in agreement with Matsuo et al [5] . The switch 
from first to second order occurs at T/T c0 — 0.195 to be 
compared with the exact result 0.076. When we go to 
iV = 3we find 0.063 and for N — 5 our results coincide 
with those of Matsuo et al at 0.076. 



Let us now consider the various steps in extending 
this method to full generality [11]. First we want to 
include additional Fourier components of the order pa- 
rameter. This is done for example by solving first the 
problem for the lowest Fourier component, then insert- 
ing the result in the self-consistency equation to obtain 
a new order parameter and iterating this process until it 
is converged. This process can be carried out for a given 
number of Fourier components, by projecting out the un- 
wanted ones. Naturally for a fixed number of Fourier 
components in the order parameter one has to general- 
ize the recursion relations Eq. (4) , but they still have the 
property that, for the physical solution, the Fourier com- 
ponents have a factorial decrease for large order. This 
makes again possible to set them to zero beyond some 
order, and to solve in this way for the other components. 
The system of equations one has to deal with is slightly 
less convenient than Eq.(4), but this is a linear system 
which is easily dealt with by standard numerical meth- 
ods. In practice, for our specific problem, we have found 
that, at the transition from normal state to FFLO phases, 
the coefficient Ai of the next harmonic cos(3qx) is of or- 
der 10~ 2 compared to the first one and produces a very 
small correction. This is in agreement with our findings 
[8] near the tricritical point. We have found this result 
for order parameters having either a single cosine or two 
cosines. We assumed it for order parameters with three 
and four cosines that we have considered. 

Next we can take the order parameter to be any lin- 
ear combination of cosines, with wavevectors in differ- 
ent directions. More generally we could take a three- 
dimensional Fourier expansion for the order parameter. 
Naturally in practice this becomes too heavy for a large 
number of components. This generalization amounts to 
replace our above index n in the Fourier expansion Eq.(3) 
by a set of integers, for example (ni,n 2 ,n 3 ) for the or- 
der parameter A(r) = 2Aicos(qi.r) + 2A2Cos(q2.r) + 
2 A3 cos(q3.r). This means that our Fourier components 
are now indexed by a set of points in a three-dimensional 
space. We make now our cut-off for \m \ + \n 2 \ + |tz3 | > N . 
Naturally we have much more components but our lin- 
ear algebraic treatment goes essentially in the same way 
and it stays still quite manageable. We note in particular 
that the values of the q^'s do not need to have any simple 
relation (in particular they are not necessarily orthogo- 
nal). 

An important point, which helps indeed very much the 
numerical calculation, is that the Fourier expansion con- 
verges quite rapidly as we already mentionned for the 
simple cosine. We have found that, for the case of two 
cosines, N = 3 was already quite enough and that N = h 
did not give any sizeable change. Actually, except for 
a single cosine, the results for N = 1 are already quite 
reasonable quantitatively. 

In this way we have been able to make some explo- 
ration of the evolution of the critical temperature when 
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the relative weights Aj of the cosines and the correspond- 
ing wavevectors vary. For example, in the case of 
two cosines, we have kept the wavevectors orthogonal or 
with an angle 7r/3 but minimized independently the two 
weights. We have found that, in cither case, the optimum 
is for equal weights. Conversely we have varied continu- 
ously the angle between the two directions, keeping the 
weights equal. Similarly for three cosines, taking equal 
weights, we have kept the directions in rhomboedral sym- 
metry, but allowed the rhomboedral angle to vary. These 
investigations make quite likely that, both for the cases 
of two and of three cosines, the stable phase is obtained 
for equal weights Aj and for orthogonal wavevectors qj 
with same lengths. In the following we present results 
only for these simple cases. 
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FIG. 1. Temperature T of the first order transition com- 
pared to the critical temperature Tfflo for the standard 
FFLO second order transition, as a function of the ratio T/p,. 
This is shown for various order parameters with one, two or 
three cosines oscillating in orthogonal directions. Actually we 
give on the horizontal axis the ratio T/T c , which is an increas- 
ing function of T/p, and is more physical. Note that, for low 
temperature, it is essentially equivalent to consider that this 
is P/Pfflo which is given on the vertical axis. 

The comparison between our results for order parame- 
ters with one, two and three cosines can be seen in Fig.l. 
Naturally the phase which is physically present is the 
one which has the highest critical temperature T /Tfflo, 
where Tfflo is the critical temperature for the standard 
FFLO second order transition. For T/T c > 0.154 the 
phase with a single cosine is the stable one, in agreement 
with previous results near the tricritical point [7,8]. On 
the other hand in the range 0.080 < T/T c < 0.154, the 
stable order parameter has two cosines. Finally the phase 
with the three cosines is the stable one below 0.080T C . 
As a result the transition from the normal state to the 
stable FFLO phase is always a first order transition in 
three dimensions, in contrast with the original proposal 
[2] of LO. We note that the succession of transitions we 
have found is in agreement with the scenario proposed by 



Houzet et al [7]. It is by no means obvious that increas- 
ing the number of cosines does not lead to even further 
increase in stability. We have looked in this direction by 
considering an order parameter with four cosines, with 
symmetrical directions pointing toward the corners of a 
cube (the angle between any two directions is 70.5°; note 
that this is an aperiodic order parameter). The result 
for T /Tfflo (1.016 at T = 0) is worst than for three 
cosines or even for two cosines. Qualitatively this goes 
in the direction proposed recently [13] of an effective re- 
pulsion between the directions of the various wavevectors 
involved in the order parameter, with an angle of repul- 
sion of order 2arccos(l/g) = 67° for q = 1.2 . Taking 
into account that A is not zero pushes this angle some- 
what above 67° and this agrees qualitatively with the fact 
that we can fit three, but not four directions on the unit 
sphere. 

It is remarkable and unexplained that T /Tfflo stays 
so near 1 for all our results, whatever the choice of the 
order parameter. However since in all cases we find that 
at the transition A /p is of order 0.2, the order param- 
eter stays rather small. This seems to imply that the 
transition is not far from being second order, which is 
consistent with the fact that T /Tfflo stays near 1. It 
would be interesting to explore if this feature is a mere 
coincidence by going to more complicated situations, for 
example by taking into account Fermi liquid effects [3] 
within the framework of Landau's theory. 
* Laboratoire associe au Centre National de la Recherche 
Scientifique et aux Universites Paris 6 et Paris 7. 
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